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Routine satellite operations for the Earth Science Constellation (ESC) include 
collision risk assessment between members of the constellation and other orbiting 
space objects. One component of the risk assessment process is computing the 
collision probability between two space objects. The collision probability is 
computed using Monte Carlo techniques as well as by numerically integrating 
relative state probability density functions. Each algorithm takes as inputs state 
vector and state vector uncertainty information for both objects. The state vector 
uncertainty information is expressed in terms of a covariance matrix. The collision 
probability computation is only as good as the inputs. Therefore, to obtain a 
collision calculation that is a useful decision-making metric, realistic covariance 
matrices must be used as inputs to the calculation. This paper describes the process 
used by the NASA/Goddard Space Flight Center’s Earth Science Mission 
Operations Project to generate realistic covariance predictions for three of the 
Earth Science Constellation satellites: Aqua, Aura and Terra. 


Nomenclature 

P(t 0 ) = constructed epoch covariance 

P(t) = predicted covariance at time t 

<D (tJo) = state transition matrix from time t 0 to time t 

Q(t) = process noise matrix 

Q mc = R1C position and velocity process noise submatrix 

Q C{ = atmospheric drag process noise matrix element 

crf lf = measured prediction error of the / th vector component in R1C frame 

A/ _ / Jlf , _A/ , _A/ 

<J ~ y&R + or [ + <J C 

erf = root-variance of the / th diagonal entry of P(t) 

. _ _ 7T 

<T r = yl&R +<?I +<?c 

I = 3x3 identity matrix 

<Tv = standard deviation of the / th component of the acceleration error vector 
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I. Introduction 

T HE Earth Science Mission Operations (ESMO) Project at NASA/Goddard Space Flight Center (GSFC) 
has established a routine conjunction assessment process for satellites in the Earth Science 
Constellation (ESC) 1 . The conjunction assessment (CA) process consists of: 

1. Predicting close approaches between members of the ESC and other objects in the United 
States Strategic Command’s (USSTRATCOM) Space Object Catalog 

2. Assessing the collision risk for any close approach predictions 

3. Performing any necessary risk-mitigating action 

Several times a week, close approach predictions are made by Cheyenne Mountain’s 1 st Space Control 
Squadron (SPCS) personnel. Close approach information is sent to NASA/GSFC where the GSFC CA 
team performs collision risk assessment analysis. One component of the risk assessment process is to 
compute the collision probability between objects. For the ESC, the collision probability is computed using 
two methods: (1) numerical integration of a relative state probability density function and (2) Monte Carlo 
simulations. In both cases, inputs to the probability calculation consist of state vector and state vector 
uncertainty information, as well as a user-specified keep-out region. 

The state information provided by 1 st SPCS is generated using observations collected from the US 
Space Surveillance Network (SSN). The orbit determination methodology employed by Cheyenne 
Mountain is a weighted batch least-squares (BLS) estimation method, which provides a “best-fit” trajectory 
(one that best matches the observations). The covariance derived through the BLS process is often found to 
be optimistic 2 * * 4 , since the errors in the dynamic model are usually not included in the covariance 
calculation. 1 st SPCS makes use of a consider parameter technique 5 , where the uncertainties of unmodeled 
systematic errors (i.e. errors that are not estimated through the orbit determination process) are still 
“considered” when the state covariance is computed. This method has the effect of adding additional state 
errors directly into the best-fit derived covariance matrix. The consider parameter method helps shape the 
epoch covariance, but the predictive covariance uncertainties are still typically underestimated when the 
covariance is propagated forward in time. 

Current operational procedures require that 1 st SPCS provide state and covariance information for close 
approach predictions that violate a specified alert threshold. After this information is received at GSFC, the 
probability calculations can be made. In addition to using the SSN-based solutions, risk assessment and 
risk mitigation analyses are performed using GSFC-derived state and covariance as inputs to collision 
calculation. In particular, this is done for the three ESC members that comprise the Earth Observing 
System (EOS); namely, Aqua, Aura, and Terra. Using mission-derived data in the risk assessment process 
allows for an independent corroboration and assessment of the uncertainty information provided by 1 st 
SPCS. Moreover, realistic covariance predictions with orbit adjust maneuvers modeled can be used to 
probabilistically assess how each maneuver will affect a given close approach prediction. 

This paper describes the process used to generate realistic covariance predictions for Aqua, Aura, and 
Terra. The process is relatively straightforward; the definitive covariance is propagated ahead in time using 
the state transition matrix and process noise matrix 0 as follows: 


P{t)=®(t,t 0 )p(t 0 )®(t,t 0 ) T +Q(t) (1) 

To ensure that the covariance predictions accurately reflect the orbital uncertainties for a given prediction 
period, calibration of the prediction errors, which is achieved by adjusting the process noise variance 
acceleration terms, must be performed. The following sections present the details of this covariance 
generation process, which consists of three steps: (1) Construction of a realistic definitive covariance 
matrix; (2) Characterization of the orbit prediction errors; and (3) Calibration of the covariance growth rate. 
Numerical results for Aura are included. 
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II. The Covariance Generation Process 

This section describes the process employed by the ESMO project to construct a realistic covariance 
prediction for member satellites of the ESC. The definitive covariance matrix is propagated ahead in time 
using the state transition matrix and process noise matrix. If the predictive covariance is to accurately 
represent the future state uncertainty, prediction errors that arise from force modeling errors must be 
included. Therefore, the definitive covariance matrices and the process noise matrix are empirically 
calibrated to match measured position prediction errors. The following sections describe this three-step 
process in more detail. 

A. Definitive Covariance Construction 

Orbit determination for the EOS missions is performed in the GSFC Flight Dynamics Facility (FDF). 
The state estimation software is the Real-Time Orbit Determination (RTOD®) system, which uses an 
extended Kalman filter to process tracking measurements from NASA’s space and ground networks. 
RTOD® is updated daily with the latest solar flux and geomagnetic index (GMI) indices. Examination of 
several months of definitive covariance solutions shows that the solutions are very consistent from week to 
week, i.e. the estimated position uncertainty numbers vary by only a few meters. Therefore, in the absence 
of any orbit adjust maneuvers, it is sufficient to create a single “representative” definitive epoch covariance 
to be used to initialize the predictive covariance generation. The representative epoch covariance P(t 0 ) is 
constructed by averaging the definitive RTOD® covariance elements using the procedure described below. 
It was observed that the correlation terms of the definitive covariance vary significantly with time - e.g. 
some of the correlations are ‘broken’ during tracking data contacts. In particular, the correlation terms for 
the in-track position/radial velocity and the radial position/in-track velocity vary considerably with time. 
The averaging method removes these short-term variations. 

Let Py(t) denote the (ij) element of the definitive covariance at instant /. Additionally, define P tJ (/ 0 ) 
as an “averaged” time-history value based on /L (/) for some specified time span. The following steps are 
performed to construct a realistic epoch covariance matrix: 

Step 1: Select a time span and compute the average value for each position/velocity variance element P u . 



Step 2\ Similar to the time series averaging performed in Step /, compute the average correlation value for 
each of the off-diagonal terms. For each (ij) entry, denote the average correlation coefficient as p fJ (f 0 ) . 

Step 3: Using the average values computed in steps 1 and 2, compute the off-diagonal terms Py{t 0 ) as 
follows: 


P,M)=Pij ('o)k('o> (3) 

For convenience, /L (/ 0 ) is abbreviated as P(t 0 ) in the remaining discussion. 

Steps 1-3 describe the construction of the 6x6 epoch covariance. The final size of is a 7x7 matrix, 

where the uncertainties due to atmospheric drag are included, since atmospheric drag modeling error is the 
primary dynamic error contributing to the prediction error of the EOS spacecraft . The drag uncertainty is 
represented by the uncertainty in a unitless coefficient Cj. Thus the 6x6 covariance is augmented to include 
an initial Q uncertainty by placing the initial Cj variance along the diagonal and the 7th row and 7th 
column off-diagonal terms are set to zero. That is, P 7 .(/ 0 )= P /7 (f 0 )= 0 for all i ^ j. 
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B. Prediction Error Characterization 

In the absolute sense, state errors cannot be measured. Sensor measurements and dynamic modeling 
errors present in the orbit determination process lead to prediction errors when the epoch state is used in 
future predictions. Thus orbit prediction errors can only be approximated by examining and trending past 
performance. The prediction error characterization methodology employed for the EOS missions is to 
measure the orbit prediction errors by comparing predictive ephemerides with definitive “truth” solutions 
over the same time period. These differences or “overlap” comparisons are used to assess the accuracy of 
state predictions. 

The truth trajectory is constructed by ‘post-processing’ the previous day’s tracking data using definitive 
solar flux and GMI indices files. Each day all measurements received up to the current time (i.e. 12:00 
GMT on the current day) are processed to create the definitive solution for the current day. Any additional 
tracking data that may have significantly lagged the filter’s real-time processing is included as well. Each 
day’s definitive state corresponding to 12:00 GMT is used as the basis for the predictive ephemeris 
generation. 

Prediction error characterization consists of performing numerous overlap comparisons between 
predictive and definitive ephemerides over some period of time (i.e. for a specified sample size). Both 
position and velocity differences are calculated and are expressed in the radial, in-track, cross-track (RIC) 
coordinate frame. Error statistics for each RIC component, including the mean and standard deviation, are 
computed every 6 hours for prediction times up to 7 days. Additional statistics, such as the root-mean- 
square (RMS) are computed as well. The RMS, which is defined as 

RMS 2 = mean 2 + standard deviation 2 (4) 

is chosen to represent the 1 -sigma measured errors. After the statistics are collected, error prediction 
polynomials are computed by fitting the RMS values in a least squares sense. In particular, a second order 
polynomial for each RIC component is computed, which represents the 1 -sigma error value for that 
component as a function of time. Each component is written as 

erf 1 - a it 2 + b t t + Cj (5) 

where ‘kf denotes the measured prediction errors for the / th component (position or velocity). 

Since prediction errors are driven primarily by unpredicted variations in the solar activity and 
subsequent errors in the atmospheric drag modeling, the polynomial coefficients must be 

routinely re-computed. This is particularly important for the in-track uncertainty since errors in the 
atmospheric drag modeling manifest themselves as in-track position errors. Thus the in-track prediction 
errors will tend to decrease as the solar flux decreases. For the other position components, radial and cross- 
track errors are primarily due to modeling errors in the Earth’s geopotential and solar radiation pressure 
(SRP). These errors tend to be more of a constant bias and typically don’t exhibit a strong ‘seasonal’ 
dependence as does drag. 

C. Predictive Covariance Calibration 

Calibration of the covariance growth rate is achieved by adjusting the process noise matrix Q until the 
propagated root-variance, denoted by erf , matches the measured sigma values erf 1 over the specified 
prediction span. The epoch covariance P(/ 0 ) is propagated forward in time using the equation expressed in 
the mean of J2000.0 (MJ2000) inertial reference frame 

p(i)= M< a)p(<oM‘aY + Q(t) ( 6 ) 

The state transition matrix <J>(^ 0 ) links the initial state at time to to the state at time t. O is a matrix of 
partial derivatives of the state at the current time with respect to the state at the initial time. The 4> matrix is 
computed numerically using the central differencing method, which is accumulated using a very small step 
size using a high fidelity model such as that listed in Table 2. 
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The process noise matrix Q(t) increases the state uncertainty at each step in an additive fashion. The 
presence of the Q matrix accounts for errors in the dynamic model, such as mismodeling of atmospheric 
drag, solar radiation pressure, and the Earth’s geopotential. The matrix Q must have the same dimensions 
as P. Q is expressed in the following block diagonal form: 


Q{t)= 


Qric 

K ^6.vl 


0l.r6 > 

Qc d j 


(7) 


where 

Q R1C — the process noise submatrix for the position and velocity vectors (in RIC frame) 

Q c = the process noise matrix element for the atmospheric drag parameter Q 


The RIC process noise matrix is a 6x6 while the Cd 4 Q’ is a singleton. The QRIC matrix is formulated 
such that the random acceleration uncertainties are directly integrated into the covariance prediction 
process. The RIC acceleration sigma values are denoted by Then the RIC root-variances are 

given as: 
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where AT denotes the propagation step-size. For convenience, Qric is expressed in the compact form 
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( 10 ) 


Values for q acc and Q C{ are adjusted so that the root-variances of of P(t) are in good agreement with 

the corresponding prediction error profiles a 4/ . In particular, comparisons between the RIC position 
differences are made, where the elements of the process noise matrix are adjusted so that 
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( 11 ) 


I erf 4 - erf I < K f for each / € {/?, /, c } and some constant K. 

The steps for adjusting the values q acc and Q C/ consists of: 

1 . choosing a set of initial values 

2. generating the predictive covariance P(t) 

3. making the comparisons described in Equation 1 1 

4. adjusting q acc and Q Cj until the ‘convergence criteria 5 are satisfied 

Changes made to the 0 C{ affect the in-track error predictions. Since errors in atmospheric drag 
modeling lead to errors in the in-track position, the Q C{ value will affect the predicted in-track position 

variance. The position variances are also directly correlated with the RIC acceleration variances; therefore 
changes to the {cr^,<xy,cF ( ; } values directly affect their RIC position counterparts. Additionally, there 

exists some correlation between the in-track prediction error and the acceleration error in the radial 
direction. Thus adjustments made to a ^ will have some effect on the in-track predictions as well. Table 1 

summarizes how elements of Q(t) map to changes in the predicted position errors when the covariance is 
propagated forward. 


Tablet: Process Noise Effects 


Q(t) Element 

Affected P(t) Position Error 

a R 

Radial, In-Track 

°7 

In-Track 


Cross-Track 

Qc, 

In-Track 


Although both Q C} and <j } - will affect the in-track error prediction, only one of them is adjusted when 
calibrating the predicted covariance in accordance with Equation 11. 


III. Aura Case Study 

The methodology described in Section II is now applied to a specific case in which a realistic 
covariance prediction is made for the Aura spacecraft. First, an ‘averaged’ epoch covariance P(t 0 ) is 
constructed. Second, the position prediction errors , af f , <Jq ] are measured for prediction times up to 

7 days. Third, the covariance prediction is calibrated by adjusting the process noise matrix Q(t) until the 
predicted position root-variances , erf , <r£ j match the measured prediction errors. 


D. Force Modeling and Initial Conditions 

Both the state and covariance are propagated forward using FreeFlyer®. Table 2 lists the perturbative 
forces modeled and Table 3 lists the state vector initial conditions used in this case study. 


Table 2: Aura Force Modeling 


Parameter 


Modeling 

Earth Geopotential Model 

30x30 Joint Gravity Model (JGM) -2 

Non-Central Bodies 

Sun, Moon (JPL DE200 Ephemeris File) 

Atmospheric Density Model 

Jacchia Roberts 

Solar Radiation Pressure Model 

Thin plate model with reflectivity coefficient 

C r = 1.4 

Effective Solar Radiation Pressure Area 

47.95 m 2 

Effective Drag Area 

47.95 m 2 

Numerical Integrator 

Runge-Kutta 8(9) with 60-second step size 
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Table 3: Aura Initial State Conditions 


Parameter 

Value 

State Epoch (MJ2000) 

16-Mar, 2006 13:19:20.000 

X(m) 

435179.356 

Y (m) 

-923269.985 

Z (m) 

6997996.863 

Xdot (m/s) 

-7079.883 

Ydot (m/s) 

-2489.657 

Zdot (m/s) 

111.781 



Mass (kg) 

3046.5 

Cj 

1.7 

Cj variance 

0.01 


E. Definitive Covariance Construction 

The epoch covariance P{t 0 ) is constructed using the averaging method described in Section II-A, where 
the position/velocity covariances are determined using definitive solutions spanning the time March 16, 
2006 - March 24, 2006. Figure 1 and Figure 2 show the definitive RIC and total position and velocity root- 
variances over the 7-day period. These figures show that there are small variations during the week, where 
the position uncertainty magnitude ranges between 1 and 4 m and the magnitude of the definitive velocity 
uncertainty values range between 1 and 4 mm/s. 






Figure 1: Definitive Position Root-Variance 
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Figure 2: Definitive Velocity Root-Variance 

After the ‘averaged’ 6x6 position and velocity covariance is constructed, an initial Cj variance of 0.01 is 
assumed. Then the final form of P(t 0 ) is given as follows, where units of the 6x6 position/velocity 
covariance are m 2 and (mis) 2 and the C d sigma is unitless: 


0.4968374877 
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F. C. Prediction Error Characterization 

Next the measured position prediction errors are characterized by applying the method described in 
Section II-B. First, for each RIC component, the RMS statistic is calculated for prediction times up to 
7 days. A least-squares polynomial fit to the RMS statistics is performed. Figure 3 shows the prediction 
error information for Aura. The blue squares on the left hand side are the measured RMS error statistics 

collected at 6-hour increments. The red line is the prediction error ‘best-fit’ polynomial erf' 1 . On the right 
hand side are the fit residuals, which are denoted by ‘plus’ signs. The residuals are the differences between 
the measured RMS prediction error and the polynomial fit. 
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□ Measured RMS Prediction Error 


Best Fit Polynomial 



Prediction Duration (days) 



+ Best Fit Residuals 



Figure 3: Measured Position Prediction Errors 

The prediction error statistics were generated using data spanning January 01, 2006 - March 10, 2006. 
Figure 3 shows that the dominant error is in the in-track direction, where the measured error grows to 
approximately 150 m after 3 days and approaches 600 m after 7 days. The measured cross-track error 
grows from 5 to 25 meters over the prediction span; the radial error is equal to about one half of the cross- 
track error. Table 4 below lists the prediction error polynomials that are plotted in Figure 3. 


Table 4: Prediction Error Polynomials 


Position Component 

Measured Prediction Error 

Radial (meters) 

<t" = -0.222 1 2 + 2.444 / + 1 .668 

In-Track (meters) 

of = 6.2 1 7/ 2 + 29.85 1/ + 5.747 

Cross-Track (meters) 

ctq = -0.203/ 2 + 4.589/ + 0.246 


D. Predictive Covariance Calibration - Aura 

After the prediction errors are characterized, comparisons are made between the measured error 
^7// , <j‘) ! , erf 7 } and the root-variance of the covariance prediction {erf ,<xf,crf }. The predictive covariance 

is ‘calibrated’ so that the differences between cr? 1 and erf are sufficiently small for each position 
component. 

To illustrate the need for the addition of the process noise matrix Q , Figure 4 shows the results of a 
comparison between the measured position error and the prediction root-variance without the addition of 

the process noise matrix. That is, P(t 0 ) is propagated forward using P(t)= d>(/,/ 0 )p(/ 0 )o(/,/ 0 ) 7 . 
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Figure 4: RIC Position Error Predictions without Process Noise 

In each of the four graphs, the lower of the two error curves is the prediction root-variance that was 
obtained by propagating P(t 0 ) forward in time without process noise. Examination of Figure 4 clearly 

shows that each root- variance |cr£ , a 1 / , ] of P(t) underestimates the measured error 9 af f } as P 
is propagated forward. There is virtually no error growth in either the radial or the cross-track components. 
Additionally, the difference between the two in-track error predictions grows to approximately 50 m after 
3 days and nearly 100 m after 5 days. Therefore, in the absence of process noise, the predicted covariance 
matrix provides a conservative or ‘optimistic’ representation of the position uncertainty. 

As described in Section II-C, the process noise calibration methodology consists of: examining the root- 
variance of the predicted covariance; comparing the predicted covariance to the measured prediction error; 
and making adjustments to Q as necessary so that the two predictions are in good agreement. Since values 
of <5 j and Q C/ both affect the in-track prediction error, it is sufficient to choose to vary only one when 

adjusting the in-track error growth. So Q C/ is set to zero and the prediction convergence is chosen such 

that, over the entire prediction span, the radial and cross-track differences are within 5 m and the in-track 
differences are within 10 m. That is, for all prediction times t, 

I _A/ P I . c I A 1 P I i a I A7 P I , r / 1 

(X ft 0 (X j <X 7 <10 CT £ cr r < 5 (12) 

Equation 12 is satisfied with the following RIC acceleration variances: 
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r \Ae- 9^ 

2Ae - 1 1 


(13) 


tfacc 


5e-9 


Figure 5 and Figure 6 show the results of adding the process noise 0-matrix to the predictive covariance 
generation. Each plot in Figure 5 contains two curves; namely erf 1 and af for each position component 
and the root-sum-square of the position components cr A/ andcr /J . 





Figure 5: Position Error Predictions with Process Noise 

The differences between the two error predictions are shown in Figure 6. The radial and cross-track 
differences remain within the specified 5 m convergence threshold; the in-track difference stays within 10 
m. This gives a total position error difference of 10 m for the entire prediction time t. 
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Figure 6: Position Prediction Error Differences 


IV. Conclusions and Future Work 

The collision risk assessment process requires realistic state uncertainty information to make accurate 
collision probability calculations. A three-step methodology describing realistic covariance predictions has 
been introduced. This method is straightforward and has been used to construct covariance prediction 
models for the Aqua, Aura, and Terra missions. In the future, this method will be enhanced to include the 
generation of realistic covariance predictions through periods containing spacecraft maneuvers. 
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